Neurolmage: Clinical 5 (2014) 332-340 




Contents lists available at ScienceDirect 



Neurolmage: Clinical 



journal homepage: www.elsevier.com/locate/ynicl 



Using longitudinal metamorphosis to examine ischemic stroke lesion /fls 
dynamics on perfusion-weighted images and in relation to final outcome 
on T2-W images 

Islem Rekik^ '^ *, Stephanie Allassonniere^ Trevor K. Carpenter^, Joanna M. Wardlaw^ 

'Division of Neuroimaging Sciences, Brain Research Imaging Centre, University of Edinburgh, UK 
''CMAP, Ecole Polytechnique, Route de Saclay, 91 128 Palaiseau, France 



CrossMark 



ARTICLE INFO 



ABSTRACT 



Article history: 

Received 22 April 2014 

Received in revised form 25 June 2014 

Accepted 21 July 2014 

Available online 1 August 2014 

Keywords: 

Longitudinal metamorphosis 
Acute/subacute ischemic stroke 
Dynamic evolution 
Perfusion 



We extend the image-to-image metamorphosis into constrained longitudinal metamorphosis. We apply it to 
estimate an evolution scenario, in patients with acute ischemic stroke, of both scattered and solitary ischemic le- 
sions visible on serial MR perfusion weighted imaging from acute to subacute stages. We then estimate a patient- 
specific residual map that enables us to capture the most relevant shape and intensity changes, continuously, as 
the lesion evolves from acute through subacute to chronic timepoints until merging into the final image. We de- 
tect areas with high residuals (i.e., high dynamics) and identify areas that became part of the final T2-w lesion 
obtained at >1 month after stroke. This allows the investigation of the dynamic influence of perfusion values 
on the final lesion outcome as seen on T2-w imaging. The model provides detailed insights into stroke lesion dy- 
namic evolution in space and time that will help identify factors that determine final outcome and identify targets 
for interventions to improve outcome. 

© 2014 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license 

(http://creativecommons.Org/licenses/by-nc-nd/3.0/). 



1. Introduction 

Stroke is the third leading cause of death in industrialized countries 
and a major cause of death worldwide. Ischemia - the commonest type 
of stroke - results from disruption of blood flow within the brain caused 
by occlusion of an artery. Efforts to improve the treatment of ischemic 
stroke target 'brain tissue at risk of infarction' (the penumbra) that 
can revert to a normal state following successful recanalization of the 
occluded arteiy within the first few hours after onset (Wardlaw, 2010). 
Multi-modality brain imaging is widely used in acute stroke manage- 
ment for diagnosis, prognosis and treatment planning. In particular, 
diffusion-weighted imaging (DWl) and perfusion-weighted imaging 
(PWI) are frequently used to detect early ischemic changes and distin- 
guish between permanently damaged and salvageable tissues. However, 
while medical imaging has improved our understanding of stroke for the 
last decades, it is still not possible to determine the 'best' treatment in 
any one individual patient (Wardlaw et al., 2012). One reason for that 
is the wide spectrum of stroke lesion patterns and the lack of under- 
standing of factors governing lesion evolution. Particularly, stroke lesions 
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are frequently multicomponent, appear and disappear, distort adjacent 
tissues and have ill-defined boundaries (Rekik et al., 2012). A second rea- 
son is the wide use of basic volumetric analysis and thresholding tools for 
assessing a stroke lesion and selecting patients for clinical trials (Davis 
et al., 2008; Furian et al., 2006; Hacke et al., 2005). These basic quantita- 
tive techniques fundamentally limit our scope for examining stroke dy- 
namics as they do not capture the direction, locality or the speed of 
stroke lesion evolution. 

This suggests the need for better methods to overcome the limita- 
tions of these techniques and accurately characterize stroke lesion 
dynamics. As highlighted in the review paper (Rekik et al., 2012), no 
patient-specific dynamic models had been developed to map in detail 
the kinetics of the lesion evolution and to estimate the spatiotemporal 
deformations that ischemic tissue undergoes. Recent work by (Rekik 
et al., 2013) began the exploration of advanced spatiotemporal (4D) 
models for studying the dynamics of stroke evolution in longitudinal 
data. The applied current deformation model enabled the identification 
of areas of contraction and expansion in stroke lesion surfaces, which 
allowed us to investigate correlations between perfusion and diffusion 
lesion surface evolutions (Rekik et al., 2013). Although this model 
proved to be a mathematically robust representation of the lesion sur- 
face, it could not incorporate image intensity measures in its abstract 
mathematical framework. Hence, there was no possibility of using this 
model to study the effects of DWl or PWI or other tissue parameter 
values on lesion dynamics. Furthermore, it was not suitable for lesions 
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where the number of components varied with time. However, scattered 
multicomponent lesions and variable patterns of perfusion and diffusion 
abnormality progression are very common (84% of our cohort) (Ogata 
et al., 2011). 

To partly overcome these limitations, we used a more versatile 
approach, the metamorphosis model — introduced by Trouve and 
Younes (2005), which handles both multi-component and solitary 
lesions and incorporates gray-scale intensities into the 3D image 
evolution. To examine the influence of perfusion values on final tis- 
sue outcome, we extend the work presented in Garcin and Younes 
(2005) from two-image based metamorphosis to time-series based 
longitudinal metamorphosis and introduce a new constraint on the 
velocity path to force its regularity in time. We identified one recent 
paper that addressed longitudinal image evolution using meta- 
morphosis. Hong et al. (2012) adapted the work of Niethammer 
et al. (2011) by proposing a metamorphic geodesic regression 
using appropriate averaging of independent pairwise metamor- 
phoses. In this paper, we present a different and more compact 
formulation for longitudinal metamorphosis using an ordered 

set a = {/o, /i In] of N images. Unlike Hong et al. (2012), 

our approach does not require a pairwise estimation of N -I- 1 
metamorphoses using geodesic regression. Instead, the method esti- 
mates a single metamorphosis that exactly meets all the intermediate 
images (obsei"vations) in 3, while enforcing regularity in time for the es- 
timated velocity field (Fig. 1). 

The remainder of this paper is stnjctured as follows. In Section 2, we 
will describe the image-to-image formulation of the metamorphosis 
theory and its extension to longitudinal data. In Section 3, we present 
the results of longitudinal metamorphosis to both synthetic images and 
perfusion ischemic lesion and investigate the relation of the metamor- 
phic residual to the final T2-w lesion boundary. Finally, we present a 
critical overall analysis, presenting the major findings and limitations 
of our model and revealing new avenues for exploration in ischemic 
stroke. 



2. Methods 

2.1. Image-to-image metamorphosis 

2.1 A. Prior to the metamorphosis era: the diffeomorphic era 

The serial medical imaging and the increasingly acquired datasets to 
study changes in anatomy and brain disease evolution have triggered 
the development of compelling mathematical frameworks based on 
deformations (Klein et al., 2009). A specific category of spatial de- 
formations has spanned the attention of researchers in this field: 
diffeomorphisms (smooth deformation with a smooth inverse). 
These were largely used in different registration models (Allassonniere 
et al., 2005; Beg et al., 2005; Holland and Dale, 2011; Klein et al., 2009) 
and became a part of the classical deformable template theory - especially 
after the establishment of the Large Deformation Diffeomorphic Metric 
Mapping (LDDMM) - pioneered by Dupuis et al. (1998) and Trouve 
(1995, 1998). The metamorphosis theory is built upon the LDDMM 
framework which is based on the idea of a diffeomorphic metric. 
This metric is a distance on the object space - seen as a Riemannian 
manifold - which results from the transportation of a metric on the 
group of diffeomorphisms by a group action. This defines a distance 
between two objects through the geodesic diffeomorphic path which 
connects one to the other. In the following, we will only consider objects 
which are images. The estimated diffeomorphic transformation g^ con- 
nects a source image Ig to a target image h as follows: Io° gT^ = h- 
The central idea ofLDDMM is thatgi is the mapping at time 1 ofadefor- 
mation path. What drives the evolution of this transformation from go 
to gi is the flow equation: given Vt(g(t)) a velocity vector field 



go = id 
where id is the identity map. 
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Fig. 1. Comparison between metamorphic regression scheme presented in Hong et al. (2012) (a) and our longitudinal metamorphosis (b). The red bold line in (a) represents the meta- 
morphic regression line determined by pairwise metamorphoses between the baseline image /(fo) and the images /,. Each individual metamorphosis is defined by the initial momentum of 
the geodesic shooting. In (b), the estimated longitudinal metamorphosis path Ji morphs the baseline image successively into (i, I2, and /a till merging with the final image I4. To avoid the 
velocity jumps we force the estimated velocity to be continuous in time at the observation timepoints {to. ti, ...}. 
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However, a major shortcoming of the diffeomorphic metric is that it 
does not allow variations in the topology of the deformed object. For 
this reason, it is commonly used to study changes in anatomical struc- 
tures (as connected sets remain connected and disjoint sets remain 
disjoint) rather than in brain lesions that may take dissimilar forms 
and undergo a topology change (such as the appearance of new discon- 
nected sets of ischemic regions in stroke or the merging of different sets 
into a bigger one). To overcome this shortcoming, a different mathemat- 
ical formulation introduced a slight modification in the idea of the 
diffeomorphic metric (Trouve and Younes, 2005). This led to the defini- 
tion of a new metric: the metamorphic metric. 

2.3.2. The metamorphosis theory 

Proposed by Trouve and Younes (2005), metamorphosis is an effi- 
cient non-linear registration method that aligns two images by jointly 
estimating an optimal velocity deformation vector field Vt and an opti- 
mal intensity scalar field /[. The idea of metamorphosis stems from in- 
cluding the infinitesimal variation of the gray level of the image in the 
metric used to estimate the geodesies of deformations accounting for in- 
tensity change and producing the metamorphic metric. The definition of 
this metric requires the use of appropriate spaces to which images, 
forces acting on them, and velocities driving the evolutions of the base- 
line image towards the target images belong. The mathematical frame- 
work for metamorphosis is composed of the three main components: 

• Image: an image / is an element of the square integrable set of func- 
tions M = L^(fl) considered as a Riemannian manifold (the object 
space) with it as the support of /. A curve (f i- /t, t g [0, 1 ]) on M is 
the evolution path of a base line image Ig. M is equipped with the 
usual metric on L^. 

• Action (force): an action g e C is a diffeomorphic transformation 
which acts upon the object space M. In the realm of a classical 
diffeomorphic deformation theory, the action g is associated to the ve- 
locity v through the flow equation ( Eq. ( 1 )). A curve ( t ^- gt, t e [ 0, 1 ] ) 
on G acting on / e M describes a path of deformation morphing / over 
the time intewal [0, 1]. 

• Velocity: for all f g [0, 1 ], the velocity field Vt belongs to the vector 
space V, which is the tangent space to the action group G. We 
adopt a similar construction as in that of Beg et al. (2005) for the ve- 
locity vector space V on which smoothness constraints are placed to 
ensure the existence of optimal smooth solutions in the space of 
diffeomorphisms for the flow equation (Eq. (1 )). We endowed the 
velocity vector space V with an inner product < .,. >v defined through 
a differential Cauchy-Navier type operator L (with adjoint L^) given 
by: </: g>v = <Af ; Lg>i2 = <Lhf, g>i2 where <. , .>^2 is the standard 
inner-product for square integrable vector fields on M and the 
Cauchy-Navier differential operator as presented in Beg et al. 
(2005) L = (-aV^ -I- y)la [la is the identity matrix). Thus, the re- 
quired smoothness of the deformations is specified by the norm of 
the space V of smooth velocity vector fields through L. 

As defined by Trouve and Younes (2005), a metamorphosis is a pair 
of curves ( (g( f ) , /( t) ) , t g [ 0, 1 ] ) on the product space . = C x M, with 
g(0) = id. The effective metamorphosis path ]{t) on M is defined as a 
combination of the action deformation path g(t) and the image path 
/(f) which is the residual of the deformation on M: j{t) = g[t) ■ l{t). 
All details of metamorphosis are available in Trouve and Younes (2005). 

As shown in Garcin and Younes (2005), the optimal metamorphosis 
curves {{g{t),l{t)),t'^ [0, 1]) - morphing the source image /□ into target 
one /i - are estimated by minimizing the cost functional U(/, v): 



U(;,v) = / \v\ldt + X 

0 u- .1 0 



Io-h = gt' 'o- This nicely brings us to the diffeomorphic flow equation 
(Eq. (1)) that is driven by a sufficiently smooth deformation vector 
field to estimate (Beg et al., 2005; Trouve, 1998). All the previous equa- 
tions have the fundamental property to remain geodesic equations in 
the image manifold M. Their solution is sought through a minimization 
algorithm, which provides an efficient and robust estimation of intensi- 
ty variation and diffeomorphisms — in case of large deformations 
(Trouve and Younes, 2005). 

This formulation unifies a common mathematical framework in- 
tensity change and deformation. The residual of the deformation 
term |6(t)|^(,j = |^+ V/t ■ Vt|^2 undertakes that the metamorphosis 
scheme includes the variations in the image evolution path /(f) induced 
by the deformation field Vf, thereby it can be viewed as a condensed 
form that sums up both variations in intensity and shape. 

As we need to minimize this energy for image matching, a 
discretization step in the time and space domains is required. We repro- 
duce the discretization scheme proposed by Garcin and Younes (2005) 
for the image-to-image metamorphosis where we approximate the 
term V/f ■ Vt by (/(t -I- 1, x -I- v(f, x)) - /(f, x)) as follows: 



,. /(t + e,x + ev(f,x))-/(f,x) a/(f,x) a/(f,x) ,^ ^ 

lim^ ■ ^ — - — ^, ' -I — ^j-^- v(t,x). 

c^o e df dx ^ ' 



(3) 



Next, we discretize the energy functional in the discrete time do- 
main of evolution [0, T] (with the size of a timestep At being such that 
T = n X At) and in the image space domain (a grid). We use a trilinear 
interpolation proposed by Garcin and Younes (2005) to define real 
values for x + Vt(x) : 

r,^/,^, : xefl^r(/f^i ) (X + Vt(x)) = /(f + 1 , X + v(t, x))eR. 



By applying the trilinear operator F to + i(x -I- Vt(x)), we get: 

V) = ErJo l^tlv +■ E'Jo E,.n ir(/,,,)(x + v,(x)) + v;, . v,\l. 

v,=v, andv^_, = v^ 



(4) 



We initialize the image evolution path from the source /q to the tar- 
get image /i through piece-wise trilinear interpolation (fixed boundary 
conditions for exact matching) : ;(f) = ( 1 - t)/o -I- 1/, ; with f e [0, 1]. Ul- 
timately, the minimal metamorphosis path composed of both the o 
ptimal image path (f ^- /t, t e [0, 1]) and the associated diffeomorphic 
path (f gt, t G [0, 1]) is calculated using the metamorphosis energy 
gradients (V/jU, Vv,Lf) in a standard alternating steepest gradient de- 
scent algorithm: 



2 



(5) 



(6) 



(2) 



2.2. Constrained longitudinal metamorphosis using N images 

Now we will present the generalization to a set of time dependent 
observations. We aim to estimate a metamorphosis (w.r.t. the meta- 
morphosis metric) for which we constrain the velocity vector field 
(vt)t e (0.1] to be continuous in time in particular at the observation 
timepoints. This is based on similar equations as before and the energy 
to minimize becomes: 



The term V/t ■ Vt represents the spatial variation of the moving image 
1[ in the direction Vt. Furthermore, the moving intensity scalar field /t is 
defined under the action of the diffeomorphism gt on a baseline image 



"(^•^) = Er:o l^tlv + Jj E':oE..nir(V,)(>= + Vt(x)) + V/,.v,|?. 
v,^ = Vf^ Vj_, = VjandVf_, = v, = Vj^, forte{l, I— 1} 
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where ( Vt)t e [o.i ] is continuous in time. We enforce the continuity of the 
velocity vector field so that the whole path is continuous in time. This is 
done as follows: for any obsereation timepoint t = tabs'- 



We choose a small time discretization step between the observa- 
tions for this definition of regularity in time to be valid. This constraint 
forces the estimated metamorphosis to follow a relevant path (w.r.t. 
our application) and makes it differ from concatenated paths as illus- 
trated in Fig. 2. 

In order to minimize this energy and get the longitudinal metamor- 
phosis using N images 3 = {/q, /i In], we exactly follow the steps of 

the gradient descent pipeline using Eqs. (5)-(6), except that we modify 
two steps. In the first step, we initialize image evolution path from the 
source image /q to the target image through piece-wise trilinear inter- 
polation (fixed boundaiy conditions for exact matching). However, we 
force the algorithm to exactly go through all obsei"vations (time-series 
images) by keeping the observations unaltered and only updating the 
intensity path connecting them. The second change is that we impose 
a time-continuity constraint on the estimated velocity which forces rel- 
evant final deformation maps in a way that would better capture the dy- 
namics of the lesions. 



3.2. MR imaging and pre-processing steps 

All MR images were acquired using a GE Sigma LX 1 .5-T MR scanner 
(General Electric, Milwaukee, Wisconsin) with a birdcage quadrature coil 
and a standardized protocol for acute stroke. The spin-echo echoplanar 
imaging diffusion tensor axial sequences and dynamic susceptibility con- 
trast echoplanar imaging PWI had 15 axial slices each of 6 mm thicl<ness 
with an interslice gap of 0.97 mm and an imaging matrix of 128 x 128 
encompassing a 240 x 240 mm field of view. MTT perfusion maps were 
generated using PWI data, full details of the image acquisition and pro- 
cessing protocol have been described previously in Rivers et al. (2007). 
The MTT and final follow-up T2 -weighted images (at > 1 month after 
stroke) were co-registered using a 3D affine transformation and their cor- 
responding visible lesions were manually delineated on every slice by an 
expert. Similarly the perfusion series were registered to the diffusion se- 
ries (Bastin and Ai'mitage, 2000). We used affine registration to compen- 
sate for these effects and account for patient motion, similar procedures 
have been applied elsewhere (Huang et al., 2008). Normalized mutual in- 
formation was used as a similarity measure since the contrast between 
the diffusion baseline image and other images is dependent upon the di- 
rection of diffusion weighting — likewise the perfusion series before and 
after the arrival of gadolinium contrast. 

4. Results 



3. Material 

3.1. Patient recruitment [10 patients] 

We test the metamorphosis model on 10 representative patients 
from a study of serial MR imaging in hyperacute stroke, representing a 
range of stroke severity (mean ± std dev.): (NIHSS = 12.6 ± 8.9), 
age (74 ± 9.47 years) and acute mean transit time (MTT) volume 
(1.78 ± 1.23-10^ mm^). We include patients who had PWI images at 
around 5 h, the second at around 5 ± 1 days and the third at 10.5 ± 
2.5 days after stroke and T2 -weighted image lesion at > month after 
stroke. Furthermore, we have checked that swelling in the recruited pa- 
tients did not distort DWl lesion boundary as it can mislead the result 
interpretations. 



Before forcing velocity continuity in time at observation times 



l^tobs discontinuous at t„f,., 
*^ New value for i'(„i„ 



1 ,- 1 I i ,. I 

"■oil., toba '-ob.i 



After forcing velocity continuity in time at observation times 




Updated velocity curve continuous In time. 



t~ t V 



Fig. 2. Enforcing the continuity in time of the estimated velocity field (v,) at obsei'vation 
timepoints. (Top) For a fixed voxel x in the image, we notice that the velocity curve v, is 
discontinuous at observation timepoints tot, as both the red and blue curves are not 
"glued" together. We enforce the continuity at the observation timepoints by associating 
a new value to v„fc at Ubs that is equal to the velocity value at tois then we update ttts to 
establish the equality between the three discrete points in time. This generates a new ve- 
locity curve (in green) that is continuous in time. 



4.1. Comparison between the proposed longitudinal metamorphosis and 
independent metamorphoses 

In this section, we compare the synthetic results of our method with 
respect to those of a modest concatenation of individual metamorpho- 
ses. Note that this concatenation does not make any assumption on 
the continuity of the deformation field so that it will most probably be 
discontinuous at the observation timepoints. This problem should be 
avoided by our method and therefore better fit the real behavior of 
the lesion. To perform this, we introduce a synthetic example where 
the comparison is easier and preliminary observations are made. 

We have tested the proposed longitudinal metamorphosis which is 
the solution of an energy minimization under constraint on synthetic 
images (two small spheres merging into one sphere then expanding 
into a bigger one. Fig. 3). Our approach is different from the estimation 
of successive independent metamorphoses from image at ti to image 
at tj then from image at tj to image at t^, etc., then gluing these bits all 
together to build the longitudinal metamorphosis scenario. More 
importantly, we show in Fig. 3 successive snapshots of the estimated ve- 
locity field magnitude for the proposed constrained longitudinal meta- 
morphosis that goes through three synthetic images (Fig. 3-B, first 
row) and for an image-to-image metamorphosis concatenating two in- 
dependent metamorphoses (one morphing image 1 into image 2 and a 
second one morphing image 2 into image 3) (Fig. 3-B, second row). The 
visual compaiison of both evolution scenarios shows a lack of a transition 
phase between metamorphosis 1 and metamorphosis II: the concatena- 
tion process of two independent metamorphoses create an abrupt change 
in the estimated velocity field. However, the constrained longitudinal 
metamorphosis includes an in-between transition phase between the 
pairs of evolution {image 1 — » image 2} and {image 2 — ► image 3} at 
the observation timepoint tj (outlined in red in Fig. 3). Therefore the evo- 
lution of the set of synthetic images seems more natural in the first esti- 
mated longitudinal metamorphosis scenario. 

4.2. Metamorphic longitudinal matching applied to perfusion MR images of 
stroke 

For every patient, we estimated a longitudinal metamorphosis of 
MTT lesion from acquisition timepoints ti to ts. Both velocity and inten- 
sity paths are presented in Fig. 4. We empirically set the trade-off pa- 
rameter a such as Jj = 0.001 for all patients. We chose Cauchy-Navier 
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Fig. 3. Comparison between the proposed constrained longitudinal metamorphosis and independent metamorphoses using the estimated velocity field. (A) Synthetic data composed of 
three images. (B) (First row) We display successive snapshots of the estimated velocity magnitude map for the constrained longitudinal metamorphosis going through image 1, then 
image 2 and finally image 3. The velocity maps at the obseivation timepoints are highlighted in red. (Second row) We visualize the estimated velocity magnitude map for an image- 
to-image metamorphosis concatenating two independent metamorphoses: a first one from image 1 to image 2 and a second one from image 2 to image 3. (C) We visualize the velocity 
vector field on the x-y plane and its magnitude map at the observation timepoints. The color bar codes the velocity magnitude: the orange color represents positive magnitude (corre- 
sponding with bright areas in the gray-scale velocity magnitude maps) and the purple color represents negative magnitude (corresponding with dark areas in the gray-scale velocity mag- 
nitude maps). 



differential operator L to provide greater smootliing of our data, justi- 
fied in our study since we focus on obtaining a quantified estimate of 
the dynamics of ischemic lesion spatial margins. We set a = 0.01 and 
y = 0.001 as suggested by Beg et al. (2005). 

We used the proposed longitudinal metamorphosis to exhibit and 
examine a single map (here residual map) generated for every patient. 
A significant property of the estimated residual map is that it condenses 
in one image both the magnitude of the deformation and the intensity 
variation in the lesion during its metamorphosis (see Fig. 5), thereby 
giving insights into the kinetic and photometric dynamic behavior of 
the perfusion lesion. We would also like to highlight that both residual 
and deformation maps were generated using the estimated in-between- 
observations It and V[ evolution paths while excluding the observation 



timepoints (our fixed boundaries) since the lesion evolution in space 
and time is mainly defined by these in-between estimated paths. 

For each patient, we reconstruct the normalized residual map 
(rMap) as: 



Forxen, rMap{x) = ^ |r(/,^i){'^ + VtW)-'t{x: 
t=o 



(9) 



Then we normalize it between 0 and 1. The estimated residual maps 
quantify the variation in perfusion values inside the lesion under the ac- 
tion of the estimated deformation field. Therefore, residual areas with 
highest values mark where the most relevant dynamic change in both 






Fig. 4. Longitudinal metamorphosis, (a) MTT maps at three acquisition timepoints superimposed with the manual segmentation of the lesion (in yellow), (b) Top row: screenshots of the 
estimated intensity path (t — /,) from t, to tj of MTT lesion metamorphosis; bottom row: screenshots of the estimated velocity path (t « v,) from ti to ts. Yellow arrows point to contracting 
areas and blue arrows point to expanding areas. Both of the displayed intensity and velocity maps are normalized (without a unit). 
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Fig. 5. Algorithm pipeline, (a) Tliree MTT images at three successive acquisition timepoints. (b) Extraction of the manually segmented MTT lesions (obseivations) that will be used in the 
longitudinal metamorphosis estimation, (c) Reconstruction of the metamorphosis-derived normalized residual map (rMap). (d) Overlaying the residual map with the final T2-w image, 
(e) Automatically thresholding the residual map and detecting the areas of highest change (in pink) — superimposed with final T2-w image. 



intensity and siiape takes place. To detect these areas with highest MTT 
variation from U to tj, we define an automated threshold as the mean 
value of rMap. Then, we generate a new normalized thresholded rMap 
by including all the voxels with values above this threshold. Fig. 5 visu- 
alizes the steps of the algorithm for every patient in our cohort. Finally, 
we compute the volumetric overlap (in %) between the thresholded re- 
sidual map and T2-w lesion w.r.t. j?na/ r2-w volume. This volumetric 
overlap indicates how much of the final T2-w lesion is occupied by the 
thresholded residual map. 

Using Eq. (9), the residual map includes the most relevant lesion 
changes in shape and intensity from fi, to t2 then to t^; however, it 
does not inform us about the direction in which the perfusion change 
is going (positive direction i.e., increasing values or negative direction 
i.e., decreasing values). To interpret the thresholded rMap with regard 
to the direction of perfusion value variation — while only focusing on 
the acute and late perfusion changes, we generated a signed MTT map 
defined as the difference image between MTT image at ta and MTT 
image at ti : 

Forxen, map4^{x) = MTr,^(x)-MTT,^(x). 



Then, we mark areas in the thresholded residual map with negative 
(blue) and positive (red) mapAMrr values. This generates a signed 
thresholded residual map (Fig. 6) that allows us to simultaneously 
look at perfusion areas that underwent the highest intensity and defor- 
mation changes with distinction of areas where MTT values decreased 
from ti to t3 (negative mapAMxr values) or increased from ti to ta (pos- 
itive mapAMTT values). Positive regions in the signed rMapthreshoided rep- 
resent 'extreme' areas where the perfusion abnormality has worsened. 
In the other hand, negative regions highlight areas where the blood 
flow bettered. 

43. Remark 

The signed thresholded residual maps do not exceed the boundaries of 
the manual outlines of MTT lesions at the three acquisition timepoints. 
We also would like to point out that one could use Eq. (9) without the 
l? norm to generate the thresholded residual map. However, we pre- 
ferred to restrict our analysis on only immediate change from acute (at 
tx ) to final (at ta) dead tissue and not to consider the intermediate change 
that is governed by many unknown variables and pathophysiological 
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Fig. 6. Signed thresholded residual maps in three representative patients. Top row: same axial slices of one MTT lesion (outlined in blue) at three successive acquisition timepoints for three 
patients. Bottom row: (a) the final T2-w image where the final dead tissue is manually outlined in red. (b) Signed thresholded residual map (red color for positive values and blue color for 
negative values) overlaid with final T2-w image. Both images (a) and (b) represent the same axial slice. 
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'laws'. Comparing only the first to final tissue states is a common clinical 
assessment routine. We used it to shed light on the thresholded residual 
map. 

4.4. Exploring the predictive potential ofMTT maps using metamorphic 
residual maps 

There exists a wide variation between patients in the total volu- 
metric overlap between the thresholded metamorphic residual 
map and final T2-w lesion between patients (mean = 45.4%, standard 
deviation = 28% median = 32.8%). Positive skewness (Fig. 7-box 1) 
suggests that most patients deviate from the median with a wider 
range of volumetric overlap. This shows that the thresholded residual 
map marking abnormal perfusion areas with highest change in shape 
and intensity identifies relatively large dead areas within the T2-w 
lesion. 

The most dynamic part of PWl lesion where perfusion was improv- 
ing (i.e., negatively signed thresholded residual map) overlapped with 
the final T2 volume with range from 0 to 67%, mean = 26.1%, standard 
deviation = 19.3%, and median = 17.1% in our cohort. This shows that 
the acute-subacute improvement of the hemodynamics of the abnor- 
mal perfusion area does not necessarily imply that it will certainly end 
up outside the final T2-lesion. The negative overlap was evenly split at 
the median of the data (i.e., zero skewness. Fig. 7-box 2). Only Figs. 7- 
box 2 and 6 show that the majority of areas in the rMapthreshoided are 
positively signed (red areas in Fig. 6), i.e., with worsened blood flow. 
This indicates that the thresholded residual map contains areas whose 
MTT values increased from ti to ta, thus, identifying areas that are 
more likely to shift into an irreversible state of tissue death. 

5. Discussion and conclusions 

In the present work, we proposed a trivial extension of the work of 
Trouve and Younes (2005) and Garcin and Younes (2005) for the 
image-to-image metamorphosis theory into a constrained longitudinal 
metamorphosis that passes exactly through the true observations, 
while forcing the geometric evolution to be continuous in time. We ap- 
plied it to serial perfusion image data obtained in the acute to subacute 
stages in patients with ischemic stroke. This provided a robust, sophisti- 
cated mathematical tool to extract both dynamic and intensity features 



from the perfusion lesion summarized in the intensity variation and 
shape deformation residual map. We then used a static T2-w image at 
> 1 month after stroke to determine the reversibility of the metamor- 
phic residual of the perfusion lesion. We aimed to provide a sophisticat- 
ed method, that is blind to clinical features (e.g., site of occlusion and 
stroke severity score), for evaluating the effect of variations in cerebral 
perfusion on tissue outcome. We hypothesized that the hypoperfused 
areas that underwent the most shape and intensity variation a) would 
have been exposed to the largest differences in perfusion values and 
b) the values would be consistent between tissues that behaved 
similarly. 

We demonstrated that the estimation of longitudinal metamorphic 
residual maps is a promising tool in tracking the spatiotemporal chang- 
es in lesion shape and intensity. This overcomes major limitations of the 
commonly used 2D or 3D voxel-based or volume-subtraction methods 
that do not allow the estimation of dynamic characteristics of lesion 
progression or regression (Beaulieu et al., 1999; Karonen et al., 2000; 
Kluytmans et al., 2000; Rekik et al., 2012; Wittsack, 2002). Our model 
is fully automated and does not require any manual landmark matching. 
It is also generic so it could be applied to other medical applications 
based on serial imaging. 

Using the automatically thresholded residual map, we showed that 
the amount of variation in MTT lesion shape and intensity identified a 
large portion of tissue that is irreversibility damaged. Thus, the MTT 
map showed promise for identifying dead tissue margins and tissue 
that suraved, although there was substantial variation in the individual 
perfusion values. This means that it is difficult to identify any one MTT 
value that differentiates tissue destined to die from tissue that will sur- 
vive between patients, as it is the dynamic properties of the MTT lesion 
that seem to determine tissue fate. The present model could be used to 
tailor more sophisticated models to predict tissue fate using the esti- 
mated evolution of the perfusion lesion shape and intensity. It could 
also be used to examine other factors that may influence ischemic lesion 
outcome such as blood pressure or pharmacological treatment. The 
signed residual maps showed that most of the abnormal perfusion 
tissue in which MTT values increased (with positive map^^^^^ values 
i.e., worsening of blood flow) belonged to the final T2-w boundary. 
However, some portions (negative overlap: mean = 26.1%, standard 
deviation = 1 9.3%) of the acutely ischemic tissue whose MTT values de- 
creased (i.e., better blood flow) also ended up as dead tissue. This 




Fig. 7. Boxplot illustrating the skewness in our cohort for the volumetric overlap between thresholded metamorphic residual map (rMapthreshoicieti) and final T2-w lesion (box 1 ). Box 2 
shows the distribution of the volumetric overlap between negatively signed rMapthresiioided and final T2-w lesion. The red line represents the median value of the volumetric overlap. 
The lower boundary of the blue box represents the 25th percentile and the upper boundary represents the 75th percentile. The upper (vs. lower) horizontal line denotes the 90th (vs. 
the 10th) percentile. 
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highlights the potential for the most active positive residuals to be used 
in further predictive models to identify the boundaries of the final dead 
tissue. 

One of the key clinical findings of our approach revealed the absence 
of a single threshold for rMap that fitted all patients. Patient-specific 
thresholds were automatically defined using the metamorphic residual 
maps. These thresholds involved both variation in MTT perfusion values 
and lesion shape change as the lesion evolved [ti, ts] — showing that the 
evolution of ischemic but still viable tissue is patient-specific (and not 
dataset-specific). This is consistent with a recent summaiy (Dani et al., 
2012) which did not identify a specific threshold or range of perfusion 
values that clearly discriminated tissue fate. Our method may have 
allowed us to capture some of the highly dynamic nature of perfusion 
values in individual parts of ischemic lesions for the first time in humans, 
as for example, have been shown in experimental models due to spread- 
ing depolarizations (Strong et al., 2007). 

Several longitudinal studies have previously evaluated lesion evolu- 
tion using standard thresholding and volumetric analysis techniques 
(Beaulieu et al., 1999; Karonen et al., 2000) to assess the combined 
role of perfusion or diffusion lesion in determining the degree of tissue 
survival or death. However, they did not explore the dynamic character- 
istics of lesion boundary evolution and their relation to its hemodynamics. 
A recent study (Carrera et al., 201 1 ) noted that MTT perfusion values - 
outside and also within the DWl lesion - could be used to improve the 
identification of final infarct boundary. Two different perfusion thresholds 
were distinctively identified for two different datasets. Interestingly, in 
our study, we showed that perfusion values and spatiotemporal changes 
from acute to later stages are patient-specific and not dataset-specific, 
thereby demonstrating that perfusion thresholds or thresholds that de- 
pend on perfusion values are specific to each patient and not to a partic- 
ular dataset. This accumulating evidence mitigates against the likelihood 
of finding universal and consistent perfusion values that identify at risk 
of infarction, dead and oligemic tissue in all patients and is expected to 
change the course of future clinical trials and patient selection criteria to- 
wards individualized medicine. 

Identifying the shift in tissue abnormality, from being 'reversible' to 
being 'irreversible' in both space and time, is still one of the main chal- 
lenges in stroke. The emergence of methods for dynamic modeling in 
stroke research shows potential for advancing our understanding of is- 
chemic tissue dynamics. The key to a nuanced understanding of how 
perfusion values influence the spatial extent of tissue that will ultimate- 
ly die or survive lies in refining the perfusion hypothesis (Butcher et al., 
2005). This hypothesis states that abnormal perfusion areas where 
blood flow improved between acute and subacute stages will recover 
and those where blood flow worsened have the greatest likelihood to 
be irreversible. However, while this may be generally true, we have no- 
ticed in our study that there is substantial variation in perfusion values 
in space and time that limits use of specific perfusion values in predic- 
tion of tissue fate. Exposing the 'laws' that govern the hemodynamics 
of stroke would revolutionize stroke research. This may be achievable 
by using a more sophisticated version of longitudinal metamorphosis 
(e.g., by including tissue heterogeneity or other perfusion parameters 
or site of occlusion) that we have demonstrated which is now feasible. 

Our study has some limitations. We did not account for differences 
in perfusion values between gray and white matter as ischemic lesions 
are difficult to segment into gray and white matter at present (Koga 
et al., 2005). We checked visually that the included patients did not 
have large amounts of lesion swelling as this could bias assessment of 
the final dead tissue volume. Accounting for acute swelling and late ex 
vacuo effect using more sophisticated registration algorithms would in- 
crease the accuracy of our method but validated techniques are not cur- 
rently available. Our sample size was small — we chose these 10 patients 
to represent a range of lesion morphologies and changes over time to 
illustrate that the method was feasible and determine its potential for 
displaying dynamic stroke lesion pathophysiology. It was not our inten- 
tion to provide definitive perfusion values or to examine how, for 



example perfusion values might influence diffusion values, or the im- 
pact of clinical variables. This is subject to future work in larger datasets 
with a more personalized version of the implemented metamorphosis 
method. 

6. Conclusion 

We used longitudinal metamorphosis to track the evolution of per- 
fusion abnormality in stroke using time-series imaging. This model pro- 
vides novel ways to identify the most active changes in ischemic lesion 
hemodynamics using signed residual maps and we believe will be valu- 
able in future stroke research to clarify what determines ischemic lesion 
evolution and identify new potential targets for development of new 
treatments to increase the change of better clinical recovery. The appli- 
cation of the metamorphosis model to scattered and solitary perfusion 
stroke lesions has led us to a set of obseivations: some are partly consis- 
tent with stroke literature such as the role that perfusion values play in 
determining final tissue fate, and others pave the way for a paradigm 
shift with regard to the 'universality' of perfusion values in determining 
final tissue fate. We have also presented a novel powerful analysis tool 
(metamorphic residual map) for extracting both kinetic and intensity 
features from time-series imaging. Future work includes integrating 
new clinical variables into the longitudinal metamorphosis model and 
validating our findings in larger datasets using different perfusion 
maps. Moreover, this model and algorithm are not specific to stroke le- 
sions and could be used to model other brain diseases' evolution with 
topology change such as hematoma or white matter hyperintensities 
or tumors. 
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